Decay of a superfiuid current of ultra-cold atoms in a toroidal trap 
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Using a numerical implementation of the truncated Wigner approximation, we simulate the ex- 
periment reported by Ramanathan et al. in Phys. Rev. Lett. 106, 130401 (2011), in which a 
Bose-Einstein condensate is created in a toroidal trap and set into rotation via a Gauss-Laguerre 
beam. A potential barrier is then placed in the trap to study the decay of the superflow. We find 
that the current decays via thermally activated phase slips, which can also be visualized as vortices 
crossing the barrier region in radial direction. Adopting the notion of critical velocity used in the 
experiment, we determine it to be lower than the local speed of sound at the barrier. This result is 
in agreement with the experimental findings, but in contradiction to the predictions of the Gross- 
Pitaevskii equation. This emphasizes the importance of thermal fluctuations in the experiment. 

PACS numbers: 03.75.Kk, 03.75.-b, 67.85.De, 67.85.-d 
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I. INTRODUCTION 

Superfluidity is a compelling and counter-intuitive phe- 
nomenon that has intrigued scientists for decades. The 
interplay of quantum motion of particles, quantum statis- 
tics and interactions gives rise to dissipationless flow, the 
defining property of superfluidity. This flow, however, 
will only be sustained within a certain parameter regime. 
If the system is perturbed by a sufficiently large pertur- 
bation, its dissipationless nature will break down. To un- 
derstand this breakdown in fact constitutes understand- 
ing superfluidity itself, as it entails understanding why 
excitations are suppressed in the superfiuid regime and 
what constitutes a sufficiently large perturbation that 
will destroy superfluidity. 

Experiments in superfiuid helium, Refs. [1, 2], seem to 
suggest that one such perturbation is an impurity or con- 
tainer wall that moves relative to the superfiuid with suf- 
ficiently large speed that it leads to a breakdown of super- 
fluidity. This indicates the possibility of a critical veloc- 
ity, above which dissipation develops and the superfiuid 
current decays. In a seminal study, Landau related the 
critical velocity to the elementary excitations [3] of the 
system. An excitation of energy e(p) with momentum p 
can only be created above the velocity v c = min(e(p)/[p|), 
while fulfilling both energy and momentum conservation. 
For a system with an excitation spectrum which has a 
roton minimum, such as helium, the excitation of rotons 
determines the critical velocity of superfiuid helium. For 
a weakly interacting system with a Bogoliubov excita- 
tion spectrum, the low-energy excitations are phonons 
with energy e(k) = hc\k\, where k is the wave number, 
and the above expression is equal to the speed of sound 
c. Feynman considered yet another type of excitation, in 
the situation where a superfiuid flows out of a channel 
into a reservoir and suggested that the relevant excita- 
tions were vortex-anti- vortex pairs [4]. Using energetic 
considerations, he estimated the critical velocity to be 
V c = [h/(md)] log(rf/a), where d is the channel diameter, 
m is the atomic mass, and a is the vortex core diam- 



eter. However, many questions about the phenomenon 
of superfluidity are still unanswered, regarding the di- 
mensionality of the system, temperature and boundary 
effects. 

With the technological advances in ultra-cold atom 
technology these questions can now be addressed in a 
widely tunable environment, in the flow of Bose-Einstein 
condensates (BECs), see e.g. Refs. [5-8]. The critical ve- 
locity that was found in [5] was much smaller than the 
sound velocity, while the ones that were found in [6] were 
comparable to it. Theoretical studies were reported in 
Refs. [9-12]. In [10] it was found that for a rectangu- 
lar barrier the critical velocity is the local sound velocity 
at the barrier, within a Gross-Pitaevskii equation (GPE) 
approach in one dimension (ID). In [12] the instability of 
the flow due to surface modes was explored. 

In the experiment performed at NIST [7] , a critical ve- 
locity less than the local sound speed was found when a 
barrier was raised into the superfiuid flow in a toroidal 
trap. Toroidal BECs, which have been proposed and 
investigated using a variety of methods [13-15], have 
recently been used to generate persistent currents and 
study their subsequent decay [7, 14, 15]. Potential appli- 
cations of toroidal BECs include high precision interfer- 
ometry [16] and analogs of SQUIDS in atomtronic circuits 
[17]. ' 

In this paper, we study the superfiuid properties of 
BECs in toroidal traps using a numerical implementa- 
tion of the Truncated Wigner approximation (TWA), 
Refs. [18-20]. This formalism includes the next order 
of thermal and quantum fluctuations beyond the CPE- 
approximation. We simulate the experiment in Ref. [7], 
and find that a GPE description is inconsistent with the 
experimental results. The TWA approach, however, sug- 
gests that thermal fluctuations are of visible importance, 
and further, it allows for the identification of the decay 
mechanism, which are phase slips resulting from vortices 
crossing the barrier region, as we discuss in this paper. 
The comparison to the experimental results suggests that 
the findings of Ref. [7] constitute 'post-GPE' dynamics, 
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in the sense that the inclusion of fluctuations is vital for 
its understanding. 

This paper is organized as follows: In Sect. II we de- 
scribe how the system is modeled in our formalism; in 
Sect. Ill we discuss the properties of the superfluid de- 
cay that we find; in Sect. IV we illustrate the properties 
of the phase slip mechanism; and in Sect. V we compare 
our results directly to the experimental measurements. 
In Sect. VI we discuss the temperature dependence of 
the decay, in Sect. VII we conclude. In Appendix A, 
we report our numerical method of determining the tem- 
perature of the ensemble. In Appendix B, we discuss 
different estimators of the chemical potential and in Ap- 
pendix C the dependence of the local speed of sound on 
dimensionality. 



II. MODELLING THE SYSTEM 

The semi-classical TWA method was developed in the 
field of quantum optics [18] and later formulated within 
a path-integral formalism [19]. In this method, an en- 
semble of initial conditions is generated from the Wigner 
distribution of the initial state and then propagated ac- 
cording to the classical equations of motion. Obscrvables 
are calculated in each realization and then averaged over. 
This method captures the next order of quantum and 
thermal fluctuations beyond GPE. Other TWA studies 
on ultracold atom systems have been reported on dipolar 
oscillations [21], non-adiabatic loading of a BEC into an 
optical lattice [22], dynamics of two-dimensional super- 
fluid bi- layers [23, 24], dynamical instabilities of a BEC 
in a one-dimensional lattice [25] and dynamics of spinor 
condensates [26]. 

To carry out the numerical simulations, it is conve- 
nient to discretize real space and represent the continuous 
Hamiltonian by the discrete Bose-Hubbard Hamiltonian 
[27] on a 3D square lattice of dimensions N x x N y x N z : 

k = - 3 E + + \ E - x ) 

(ij) * 

i 

where J is the hopping parameter, U is the on-site energy, 
ipj(ipj) are the bosonic creation (annihilation) operators 
at site j, and (ij) indicates nearest-neighbor bonds. For 
a lattice discretization length I, the Bose-Hubbard pa- 
rameters are related to the continuum parameters, cp. 
Ref. [28], by J = h 2 /(2ml 2 ) and U = gl~ 3 , where m 
is the atom mass, and g = 4TTa s h 2 /m; a s is the s-wave 
scattering length. The real space location r = (x, y, z) 
is related to the lattice location i = (i x ,i yi i z ) through 
r = li. L x = IN X , L y = IN V , L z = IN Z are the dimen- 
sions of the discretized space used in the simulations. To 
represent approximately the toroidal geometry we use pe- 
riodic boundary conditions along the x-direction. The y 



direction represents the radial direction and the z direc- 
tion is the direction perpendicular to the plane of the 
torus, which corresponds to the vertical direction in the 
experiment. The origin is located at (L x /2, L y /2, L z /2). 

The time-dependent external potential Vi(t) = 
Vtr.i(t) + Vb,i(t) consists of the harmonic trap, V tr ,i{t) = 
a(t) (oj 2 y 2 + uj 2 z 2 ) /4JZ 2 , with trapping frequencies uj v 
and uj z and a Gaussian barrier potential, Vbs{t) — 
v(t)Vbo e x P [— (x — xi,) 2 /2l 2 ] . Vbo is the strength of the 
barrier, Xb is its location, and h its width. The time- 
dependent coefficients, a(t),r](t) are varied between and 
1: The trapping potential is ramped up adiabatically to 
create the initial state, as described below, and the bar- 
rier potential is ramped up similarly to the experimental 
procedure. As discussed in Ref. [28], the discrete model 
approximates the continuum system when the healing 
length £ = y^h 2 /mgn ^ at some location i and the ther- 
mal de Broglie wavelength, A = y // 2i:h 2 /mksT arc com- 
parable to or larger than the lattice spacing I, where 
n o,i — (ni) is the density, ks is the Boltzmann constant, 
and T the temperature. 

Within the TWA approach, the operators ipi (t) are re- 
placed by classical fields tpi(t), which propagate accord- 
ing to the equations of motion derived from Eq. 1. We 
initialize the fields ipi(t = 0) according to the Wigner 
distribution of a homogeneous (i.e. Vi(t = 0) = 0), 
weakly interacting Bose gas, within the Bogoliubov ap- 
proximation, as in Ref. [24] . Using the Bogoliubov trans- 
formation in the phase-density representation [28], the 
Hamiltonian (1) is mapped to H' — "52, v e v $ v c v , where 
e v is the energy and c\,(c v ) are the creation (annihila- 
tion) operators for the Bogoliubov modes. Terms be- 
yond quadratic order in the fluctuations are ignored. 
The Wigner distribution for a thermal ensemble of har- 
monic oscillators at temperature To is a product of Gaus- 
sians, W <~ \[ v e - x v/ 2cr x,v-Pv/ 2 °p,v w ith variance a\ = 
l/[2e y tanh(e I/ /2T )] and a 2 p v = e v / [2 tanh(e^/2T )i for 
the position and momentum, which are mapped to the 

Bogoliubov modes using c v = (i/y/2e v ) p v +(^\Je v /2^ x v . 

After the initialization, the harmonic trapping potential 
is slowly ramped-up to generate the ensemble in the trap. 
We next measure the temperature of the trapped ensem- 
ble, as described in Appendix A. After this the experi- 
mental procedure of Ref. [7] is simulated. We imprint a 
phase winding -0(x) — > e~ l ^ x ^(x), with tp(x) = 2ttx/L x , 
then ramp up the barrier, hold it constant for approxi- 
mately 2s, and ramp it down. 

In Fig. 1 we show an example of this process to demon- 
strate the preparation sequence. In this particular exam- 
ple, we choose N x = 126, N y = 21 and N z = 7; n = 8.1 
is the expectation value of the density of the initial ho- 
mogeneous ensemble; To/J= 1 is the temperature of the 
initial state. We set U = 0.07 J, which translates into a 
lattice spacing of I — 0.987 /urn, a time step At = h/J = 
0.706 ms and energy scale J = 10.8 nK for sodium atoms. 
In Fig. 1(a) we show the homogeneous density in the 
x — y plane at z = and the azimuthal current, defined 
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FIG. 1. (Color online). Density and current in ^-direction in 
the x — y plane at z = 0. Panel (a) shows the homogeneous 
ensemble after initialization at t n — lOOft/J; (b) shows the 
ensemble after ramping up the trapping potential at i„ = 
27080ft/ J; (c) shows the system after phase imprint at t n = 
27120ft/ J and (d) after ramping up a barrier, located at x — 0, 
at t n — 28000ft/ J. The total lattice size is given by N x = 
126, N y = 21 and N z = 7. The total atom number is TV = 
150028. The barrier strength is V b0 /J = 2.65. The ensemble 
is initialized with To/ J = 1; after turning on the trap the 
temperature is T = (5.48 ± 0.11) J = (59.2 ± 1.2) nK. The 
data shown here is averaged over 512 realizations. 



the barrier strength is V&o = 2.65J. The barrier is held 
at its maximum height for 2850ft/ J « 2 s and ramped 
back down linearly over Atb ■ These time scales are based 
on the experimental procedure. The density depletion 
at the barrier is apparent in Fig. 1(d). Simultaneously, 
due to the constriction at the barrier, the current at 
z = increases at the barrier, while the total current, 
i.e. the current integrated over y and z, is unchanged. 
The total "experiment" time following the initialization 
is 3600ft/ J fa 2.5 s. All the numerical results presented in 
this paper use the same lattice discretization and times 
described here. 

In the following, we discuss numerous simulations, in 
which the parameters of the system are varied. The to- 
tal number of atoms ranges from 50,000 to 180,000. The 
number of lattice sites in y- and z-direction are chosen 
to be larger than the Thomas-Fermi radii of the conden- 
sate in these directions, and therefore vary with the total 
number of atoms. L y ranges from 17 — 21Z and L z ranges 
from 5 — 71. In the elongated direction, the length is 
L x = 12QI = 124.4 fim. For a ring with circumference 
L x , the radius is R = 19.8 /xm. 



as jx{p) = —iJlh 



at z = 0. The 



current j x (r) has a zero expectation value in the initial 
state. Both quantities are averaged over 512 realizations. 
Next, we slowly ramp up the harmonic trap, Vt r ,i{t n ), ac- 
cording to a(t n ) = {1 — tanh[(£„ — iotr)/Ttr]} /2, where 
T tr = 3200/yj = 2.26 s and t otr = 8100h/J = 5.72 s. 
The trapping frequencies in the y— and z— directions are 
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0.5J = 2tt x 113 Hz, and hw z = 2.5J = 2tt x 563 



Hz. In Fig. 1(b) we show n(x, y, 0) and j x (%, y, 0) at time 
tn = 27080/t,/ J after the trap is fully ramped on and the 
system has been allowed to equilibrate. The initialization 
in this example is completed at 27000ft/J. We introduce 
the time variable t = t n — 27000ft/ J, which corresponds 
to the time of the experiment, while t n corresponds to 
the numerical time, including the initialization process. 
After the trap has been ramped up, the density is in- 
homogeneous and has a maximum at the center in in- 
direction. Although the current has a zero expectation 
value, some fluctuations are visible due to the finite tem- 
perature that has been introduced by the initialization 
process and the trap ramp-up. Using the temperature 
measurement of Appendix A, we determine the temper- 
ature to be T — (5.48 ± 0.11) J for the example in Fig. 1. 

Next, we imprint a 2-7T phase winding at time t p h — 
27100ft/J. The density and current just after phase im- 
printing, at time t n = 27120ft/ J is depicted in Fig. 1(c). 
The density profile is unchanged, but the current now 
has a finite expectation value, as the atoms circulate to 
the right, and displays some thermal fluctuations. A bar- 
rier potential is ramped up at 300ft/ J ps 0.2 s after phase 
imprinting. The barrier is centered at Xb = 0, has a 
1/e 2 width of 2lb = 61, and is ramped up linearly as 
v(tn) = (*n — tob)/Atb over a time Atb = 145ft/ J ~ 0.1s 
starting at iofc = 27400ft/ J. For the example in Fig. 1, 



III. DECAY OF CURRENT & CRITICAL 
VELOCITY 

In the absence of the barrier potential Vb(x) we find 
that the superfluid circulates essentially without decay 
on the simulation time scales, consistent with the exper- 
imental findings. However, for non-zero barrier heights, 
decay can occur. In order to characterize this decay, we 
define the azimuthal component of the average total cur- 
ieatj T = (N 9 N y N s )- l ^ r j x (T). 

In Fig. 2 (a) the average total current, jx, is plot- 
ted as a function of time for different barrrier heights, 
for N = 51408 atoms on a 126 x 17 x 5 lattice. The 
barrier height is reported in units of the bulk chemical 
potential in the absence of the barrier, which is calcu- 
lated from the density at the trap minimum, averaged 
over the azimuthal direction, which for this example is 
/io = 2.14J (see Appendix B). The barrier is ramped on 
and off linearly during the time indicated by the gray- 
shaded regions. The temperature of the TWA simula- 
tions is T = (4.34 ± 0.08) J = (46.9 ± 0.9) nK, measured 
before the barrier ramp up, as described in Appendix A. 
For comparison, we also show results from GPE simula- 
tions in Fig. 2 (b). For the GPE simulations, the initial 
state was generated by using imaginary time propagation 
to calculate the GPE ground state in the trap [29]. In 
the GPE simulations, two scenarios are observed: Either 
the current decays quickly compared to the experimental 
time scale, or it persists, and the different barrier heights 
only effect the time-averaged value of the current while 
the barrier is up. For each case an oscillatory behavior is 
observed that has little damping on the time scale of the 
experiment. These oscillations are due to excitations gen- 
erated during the barrier ramp-up and the amplitude of 
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FIG. 2. (Color online). Time evolution of the average total 
current jr for various barrier heights. The simulations shown 
here use a 126 x 17 x 5 lattice, and N = 51408 atoms. The 
barrier is ramped on and off linearly during the times marked 
by the shaded regions. Panel (a) shows TWA simulations, 
which use 256 realizations. Prior to the barrier ramp-up the 
chemical potential is no = 2.14J = 0.86&j z and the temper- 
ature is Ttwa = (4.34 ± 0.08) J = (46.9 ± 0.9) nK. Panel (b) 
shows GPE simulations, with no = 2.2 J — 0.88hia z . 



these oscillations decreases as the barrier ramp-up time, 
A tb , is increased. In the TWA simulations the decay be- 
havior crosses over smoothly from small to large barrier 
heights. For small barriers, the decay is much slower 
than the experimental times. For larger barrier heights, 
superfhiid decay is visible, as can be seen in the exam- 
ples for Vbo/fifQ = 0.55 and Vbo/Vo = 0.59 in Fig. 2 (a). 
As the barrier height is further increased, a fast decay 
is visible, as in the example Vbo/Mo = 0.63 in Fig. 2 (a), 
that is qualitatively similar to the fast decay visible in 
the GPE simulations. We note however that this decay 
sets in at smaller values of Vbo, and thus GPE overesti- 
mates the stability of the superfhhd flow. We also note 
that within the TWA simulation, the oscillatory behavior 
of jt is damped. 

In Ref. [7] an experimentally motivated definition of 
the critical velocity of the superfluid was introduced: Af- 
ter the hold time of the barrier of typically 2s, the bar- 
rier was removed and the atomic cloud was allowed to 
expand. From the time-of-flight images it was deduced 
whether there was phase winding 1, indicating persist- 
ing superflow or phase winding 0, indicating decay of the 
superflow. For a given barrier height and chemical po- 
tential, these events occur probabilistically. The critical 
velocity was the velocity at the barrier of the initial su- 
perfluid state, which had a chemical potential such that 
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FIG. 3. (Color online). Velocity and speed of sound at the 
barrier versus barrier height in units of the bulk chemical po- 
tential. GPE data: velocity and speed of sound from the 
one-dimensional (ID) integrated density and current. TWA 
data: velocity and speed of sound from the column density 
and column current (see Appendix B). Vertical lines indicate 
the standard deviation of the velocity and sound speed. In- 
set: enlarged view, close to the GPE critical barrier height. 
The simulation parameters of the TWA data are the same 
as for Fig. 2: The atom number is N = 51408, the chem- 
ical potential no = 2.14J = 0.86hoj z , and the temperature 
Ttwa = (4.34±0.08)J = (46.9±0.9) nK. We use a 126 x 17x5 
lattice. 



half of the initial states of phase winding 1 decayed to 
phase winding 0. Phrased differently, the velocity was 
called critical if the superfluid decay time equaled the 
hold time of the experiment. A full comparison to the 
experiment will be given in Sect. V, where we imitate 
the experimental procedure. 

Here, as a first comparison, we calculate the velocity at 
the barrier maximum at the end of the hold time for the 
TWA simulations. The velocity at the barrier maximum 
is defined as 



V2D 



jx-2o(n) 



\/n 2 D(rb)n 2 D{rb + Zx) 



where r h = (x b ,0), j x -2D{*b) = J2 z 3x(xb,0, z) is the col- 
umn current and n2D{^b) = XL n ( x b, 0, z) is the column 
density at (xb, 0). We note that j x ( r b) describes the cur- 
rent along the bond (r&, rf,+Zx). In Fig. 3, this velocity is 
plotted as a function of the barrier height in units of /iq. 
For the TWA simulations, v 2 d is averaged over a time 
window t G [3250, 3395]/i/ J immediately prior to ramp- 
ing down the barrier, which is a good measure of the ve- 
locities just before the time-of-flight measurement in the 
experiment. By comparing the 'local' chemical potential, 
jj,(x b ) (see Appendix B) to the confining energies fno y and 
hw z , we find that the dynamics in the z-direction is frozen 
out at the barrier, while the Thomas-Fermi approxima- 
tion holds for the radial direction, which is why we choose 
to plot v 2 D- With these assumptions the phonon veloc- 
ity can be calculated, and expressed in terms of the peak 
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density, leading to the approximate local speed of sound 
[30,31], 

C2D(xb) = \/2^2D{x b )/im (2) 

where H2D{x b ) is the chemical potential for a quasi- two- 
dimensional(2D) system (see Appendices B and C). The 
local speed of sound is also plotted in Fig. 3. 

In the GPE simulations, the dynamics at the barrier 
in both y and z direction are frozen out, making the sys- 
tem locally quasi-one-dimensional (ID), when the barrier 
height Vbo approaches the critical barrier height. Thus 
we plot the ID velocity 

_ jx-w(n) 

\f ni D {Y b )ni D (T b + lx) 

based on the integrated ID current, jx-iD^b) = 
12 y z 3x(xb, V, z) and the integrated ID density, 
^ii)( r b) == Yly z n ( x bi V, z )- The ID speed of sound 
at the barrier is also plotted, 

cw(xb) = \/niD{x b )/m, 

where fJ-iD{x b ) is the local chemical potential for a quasi- 
1D system (see Appendices B and C). For the GPE sim- 
ulations, the velocity is averaged over a window in the 
range t g [545, 3395]ft/ J, during which the barrier is at 
its maximum value. As noted before, in the GPE simu- 
lations, no sizable decay was observed, allowing for this 
long time interval to be used for averaging. The verti- 
cal bars represent the standard deviation of the barrier 
velocity and local speed of sound, which indicates the am- 
plitude of oscillations occuring in these quantities. These 
are quite large for the GPE simulations and are due to 
the undamped oscillations generated during the barrier 
ramp- up seen in Fig. 1. 

For the TWA (GPE) simulations, the system is quasi- 
2D (ID) at the barrier only when the barrier height is 
close to the critical barrier height, but we choose to plot 
2D (ID) quantities for all barrier heights because we are 
most interested in the critical region. We note that the 
healing length at the barrier is £ « 2.2 — 2.31, which 
is shorter than the total width of the barrier 2l b = 61. 
We thus expect that the system is still consistent with a 
local-density approximation. 

We recognize in Fig. 3 the features of the TWA and 
the GPE simulations found before: As the barrier is in- 
creased, for both simulations v b increases gradually, and 
then falls off to zero as the superfluid flow becomes unsta- 
ble. The transition occurs over a finite range in the TWA 
simulations, while within GPE, there is a sharp jump be- 
tween barriers for which the current persists and those 
for which it decays. As can be seen from Fig. 3, the crit- 
ical barrier height within TWA is lower than for GPE 
and the transition from persistent to decaying current 
occurs when the local velocity is measureably lower than 
the local speed of sound, while we find (see also the in- 
set) that within GPE the current decay occurs when the 
velocity at the barrier is comparable to the local speed 
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FIG. 4. (Color online). Average total current, ]t for ten in- 
dividual TWA realizations (gray); jt averaged over 10 shown 
realizations (shown as light blue); jr averaged over 256 re- 
alizations (dark blue). The atom number is AT = 150028, 
the barrier height Vbo ~ 2.7 J = 0.68/in, and the temperature 
T = (5.58 ± 0.15) J = (60.3 ± 1.6) nK. 



of sound, as has been found in other GPE simulations in 
the hydrodynamic limit [10, 11]. 

These results suggest that within GPE, the excitations 
that lead to decay of the current are phonons, while other 
decay mechanisms are present in the TWA simulations. 



IV. PHASE SLIP DYNAMICS 

In order to understand the decay mechanism in the 
TWA simulations, we investigate the dynamics of the cur- 
rent and phase for individual realizations. In Fig. 4 the 
average total current j't is plotted for several individual 
realizations of TWA simulations (gray lines) . Within an 
individual realization there is a rapid transition from a 
state with circulation to one without circulation, while 
the average over many realizations leads to the expo- 
nential decay for the ensemble. The time when the de- 
cay occurs is probabilistic and is governed by the decay 
timescale for the ensemble. To further understand the 
mechanism of current decay, we look at the phase dy- 
namics of a single realization around the time the to- 
tal current decays. We define the phase <j>(x) along the 
density maximum, and recall that the gradient of the 
phase is related to the velocity by v = ;^V0. On a 
lattice, the phase can thus be calculated as 4>{x, 0, 0) = 

Ex'=o sin_1 where we set <M°>°,0) = 0. The 

phase along the center of the ring, <fi(x, 0,0), is plotted 
in Fig. 5 for a time window after phase imprinting and 
barrier ramp-up. Initially there is a global phase winding 
of 2ir and the steepest slope of the phase occurs across 
the barrier, where the density is at a minimum and the 
velocity at a maximum, due to flow continuity. Around 
t = 1115ft/ J the phase at the barrier jumps sharply and 
the overall phase winding drops to zero. This coincides 
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FIG. 5. (Color online). Time evolution of the phase at the 
center of the ring, <j>(x, 0, 0) versus x and t for a single TWA 
realization, with the same parameters as in Fig. 4. In Panel 
(a) we show the time interval [1025, 1225]ft/ J, in Panel (b) the 
interval [1100, 1130]ft/ J, as a contour plot. These intervals 
bracket the phase slip event at t « 1115ft/ J. 



with the total current dropping to zero. Subsequently 
long wavelength excitations are observed in the phase as 
the system dissipates the energy generated by the phase 
slip. In Fig. 5(b) large fluctuations in the total phase 
winding occur around the barrier at times close to the 
phase slip. The time for the superflow to decay in this 
example is approximately 5h/ J « 3.5 ms. 



For the same realization, several snapshots of the cur- 
rent field j X y[x, y, 0) and the out-of-equilibrium density 
5n(x, y, 0) = n(x, y, 0) — n(x, y, 0), are plotted in Fig. 6, 
where n(x, y, 0) is the time-averaged density over a time 
interval of 400fi./J. The data is averaged over five time 
steps to reduce the high frequency fluctuations. Initially, 
the current flows to the right (a). In frames (b), (c), and 
(d), a vortex in the current, accompanied by local de- 
pletion of the density, is observed in the barrier region. 
This vortex crosses the barrier region, leading to the 2ir 
phase slip observed in Fig. 5 and resulting in decay of 
circulation in the ring. After the phase slip, locally the 
current flows to the left and a density pulse moving to 
the right is observed in frames (e) and (f). Due to the 
fluctuations in the system, the vortex wanders around 
in the barrier region before crossing to the other side, 
causing the jumps in the phase at the center of the ring 
observed in Fig. 5(b). 



t=1106ft/J (d) 



t=1 114 ft/J 




FIG. 6. (Color online). Current jx, y (x,y,0) (shown as a 
vector field) and density difference Sn(x,y,0) from the av- 
erage density (colormap) in the barrier region for the same 
TWA realization shown in Fig. 5. Data plotted for a cut 
at z = for times (a) t = 11W.5H/J, (b) t = 1111.5/i/J, 
(c) t = 1112.5ft/J,(d) t = 1114ft/ J, (e)t = 1117.5ft/ J, (f) 
t = 1122.5ft/J. Each plot is averaged over 5 timesteps, which 
corresponds to a time window of 1.25ft/ J. 



V. COMPARISON TO EXPERIMENT 

In this section we compare the simulations directly to 
experiment. First, we analyze the relationship between 
the critical barrier and chemical potential, and second, 
the critical velocity compared to the speed of sound. In 
order to imitate the experimental procedure, we deter- 
mine the critical barrier height as follows. We fit the 
decay of the total current with an exponential function 
Jt = Jo cx p[~ (t — to)/r], with jo and t being fitting pa- 
rameters, to determine the decay time scale r. As dis- 
cussed before, the analysis of the experiment determined 
the critical parameters by finding the conditions under 
which half of the initial realizations had decayed to zero 
phase winding. Therefore we define the critical decay 
time scale as r cr = At/ log 2, which corresponds to 50% 
probability of decay of the current during the experiment. 
Here At = 2850H/ J is the time the barrier is held at its 
maximum value, and we ignored the decay which occurs 
while the barrier is ramped on and off. We then interpo- 
late the r- versus- Ho relation with t = tq cxp(aVho), and 
define the critical barrier height as Vf, cr = Vb(r cr ). For 
the GPE simulations, the critical barrier value shown is 
the average of the largest barrier height with no decay 
and the smallest barrier height that shows decay. 
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FIG. 7. (Color online). In Panel (a) we show the chemical po- 
tential versus the critical barrier height, based on two ways of 
analyzing the simulation results. The open symbols depict the 
chemical potential (io, calculated from the density distribution 
prior to the barrier ramp-up, and the critical barrier height 
P = Vbcr, using the barrier potential directly. The solid sym- 
bols depict the analysis that resembles the analysis of Ref. [7] . 
We use (in as the approximation for the chemical potential, 
calculated from the Thomas-Fermi distribution for N atoms, 
and p — (in ~ (iiD{xb) as the approximation of the barrier 
height, where (iioixb) is the local chemical potential at the 
barrier, determined by the column density. The temperature 
of the simulation is Ttwa = (5.50 ± 0.13) J = (59.4 ± 1.4) 
nK. In Panel (b) we show the comparison of the experimental 
data and the simulation data processed in the same way as 
the experimental data in Ref. [7]. The experimental data are 
the gray symbols, the simulation data are the solid, blue and 
purple symbols, which are the same as in Panel (a). 



In Fig. 7(a) we show the density-based approximation 
for the chemical potential (see Appendix B) plotted as 
a function of the critical barrier height Vb cr , depicted by 
the open symbols for the TWA and the GPE simulations. 
For fixed chemical potential, the current decays for lower 
barrier heights when quantum and thermal fluctuations 
are included, i.e. in the TWA simulations, compared to 
the GPE simulations. These findings are consistent with 
the results shown in Figs. 2 and 3. 

In the analysis of the experiment, the chemical po- 
tential was calculated from the total number of atoms, 
assuming a Thomas-Fermi distribution - in the absence 
of the barrier - in both radial and z-direction, (jln = 

1/2 

{gNmLjyUJ z /(irL x )] , where N is the total number of 
atoms. The barrier height was approximately deter- 
mined by j3 = (In — (i2D{xb) where (i,2D{xb) is the 
chemical potential at the barrier maximum, (i2D{%b) = 
mui z / 2nhgn2D{xb, 0), and was determined from the lo- 
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FIG. 8. (Color online). Critical velocity versus local speed 
of sound, C2D = \/2/x(x{,)/3m. The experimental data (from 
[7], depicted as gray) is rescaled to C2D- The gray shaded 
region represents the estimate of the critical velocity for tem- 
peratures in the range of 10, 20 and 30nK, as described below. 
The experimental data is inconsistent with v c /c2d = 1- Ther- 
mally activated phase slips offer a possible explanation of the 
reduction of v c /c2d below 1. 



cal column density, ri2D{xbi 0) [32]. As discussed in Ap- 
pendix B, this expression for the chemical potential as- 
sumes that at the barrier only the harmonic oscillator 
ground state is occupied in the vertical direction, which 
is valid when (i(xb) < huj Zl as is the case for all of the 
data presented here. 

In order to compare our results with the experimental 
results, we generate quantities similar to those studied 
in the experiment. In Fig. 7(a) we plot the chemical 
potential estimate (In and the approximate critical bar- 
rier /3 = /ijv — (J-2D(xb) as approximate quantities based 
on the TWA and GPE data, depicted by solid symbols, 
in comparison to the results for (1q and Vb cr - We note 
that the approximate quantities overestimate the criti- 
cal barrier potential. In Fig. 7(b) the simulation results 
for the chemical potential estimate (In and the approx- 
imate critical barrier j3 = /ijv — H2D{xb) are re-plotted, 
along with the experimental data. In these plots, we see 
that the experimental results are close to both the TWA 
and GPE simulations. As a tendency, one can speculate 
that the GPE simulations predict decay at barriers larger 
than those observed in the experiment, while the TWA 
simulations are closer to the experimental results, but the 
difference is not significant compared to the experimental 
error bars. 

However, we do find significant differences for the crit- 
ical velocity, which we present now. 

One finding of the experiment was that the critical 
velocity - as defined in Ref. [7] - is less than the speed 
of sound at the barrier maximum, which was determined 
from the column density. In Fig. 8 we show a re-analyzed 
version of the data of Ref. [7]; in particular, the critical 
velocity, normalized by the local speed of sound, is plot- 
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ted as a function of the local speed of sound [33] . 

To calculate the critical velocity for the TWA sim- 
ulations, the realizations were divided into two groups 
at each time step: those in which unit phase wind- 
ing persisted, and those in which the current had de- 
cayed. The persistence or decay of the current was 
determined from the phase winding around the ring, 

ct>(L x , 0, 0) -0(0, 0, 0) = £"=0 x sin" 1 (^) ; wh i c h was 

calculated at each time step. Subsequently the current 
at the barrier was calculated by averaging only over the 
realizations in which the global current remained. 

For the GPE simulations, we plot the mean velocity 
based on current and density integrated radially, «i_d, 
which typically differs from v^d by less than 5 percent. 
The shaded bars on the GPE data represent the mag- 
nitude of oscillation of the velocity during the dynam- 
ics. For the TWA simulations, data is presented with 
T = 5.50J = 59 nK, as well as one data point with 
T = 4.34J = 47 nK. Due to the heating that occurs in 
the initialization scheme used for the TWA simulations, 
data sets with lower temperatures were not generated. 

Interestingly, the comparison between the experimen- 
tal data and the simulations, shown in Fig. 8, demon- 
strates that the experimental data is not consistent with 
the GPE simulation. In the GPE simulation, a critical 
velocity of the magnitude of the local phonon velocity 
is predicted, however, the experimental findings are sig- 
nificantly below this value. However, the TWA simula- 
tions offer a possible explanation: Thermally activated 
phase slips, which can also be visualized as vortices pass- 
ing through the barrier, lead to a reduction of the critical 
velocity. In the experiment, the temperature is estimated 
to be of the order of 10 - 20nK, but slightly higher val- 
ues cannot be excluded [34] . The TWA simulations at 59 
nK underestimate the critical velocity of the experiment, 
while the data point as 47 nK is within the range of the 
experimental error bars. This trend is consistent with 
the expectation, that the experimental temperatures are 
lower than those used in the simulations. The shaded 
regions represent an estimate of the critical velocity for 
temperatures in the range of 10, 20 and 30nK, based on 
extrapolating the TWA results for /i = 2.14J, as we 
discuss in the next section. We also note that disorder 
of the trapping potential of Ref. [7] will also reduce the 
critical velocity found in the experiment, in addition to 
the thermal noise. 



VI. TEMPERATURE DEPENDENCE OF 
SUPERFLUID DECAY 

In this section we study the temperature dependence 
of the superfluid decay. We keep the system parameters, 
such as the barrier height, fixed and vary the temperature 
only. In Fig. 9(a) we plot the total current as a function 
of time for different temperatures, with the barrier height 
fixed to Vbo = 1.16J = 0.49/io- The total current is then 
fitted with the fitting function jx = Aexp(— t/rj), over 
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FIG. 9. (Color online), (a) Average total current jx versus 
time for different initial temperatures and numerical fits with 
Aexp(— t/r) for fixed barrier height Vbo = I.16J = 0.49/io- 
The total atom number is N = 51407 atoms, the lattice di- 
mensions are 126 x 17 x 5. (b) Time scale r of the current 
decay versus inverse temperature on a log-linear scale. Nu- 
merical fits yield Tj = (0.155±0.046) exp [(61.3 ± 1.6) /T] and 
r 3 = (6.23 ± 4.52) x lour- 10 ' 3 * ' 41 . 



the time window that the barrier was fully ramped on, to 
determine the decay time scale Tj. In Fig. 9(b), we plot 
this time scale as a function of the inverse temperature 
on a log-linear scale. We observe a strong temperature 
dependence. 

To quantify this, we attempt to fit the data with both a 
power law and an exponential function. The exponential 
scaling is motivated by theoretical work on superfluids 
[2] and thin superconducting wires [35] , which found that 
the timescale for a 2-7r decrease in the phase winding has 
an exponential dependence, r ~ exp (Af/fc^r), where 
AF is the minimum free energy barrier between the two 
states. The algebraic scaling is motivated by the behav- 
ior of ID superfluids, see e.g. [36]. We note that at the 
barrier the system is close to quasi ID, because the mean- 
field energy is smaller than uj z and comparable to uj v . As 
we see in Fig. 9(b), the data is consistent both with ex- 
ponential scaling as well as with power-law scaling. It is 
not possible to distinguish between the two because the 
temperature range that is accessible experimentally and 
in the simulations is rather narrow. We also note that 
the physical setup might not give rise to pure exponen- 
tial or algebraic behavior, because phase slips can both 
occur in the region at the barrier that is close to quasi- 
1D and also away from that region. However, the strong 
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FIG. 10. (Color online), (a) Critical barrier height normal- 
ized by the chemical potential, as a function of the initial 
temperature for TWA (circles) and GPE (square) simula- 
tions. A linear fit to the TWA data only gives Vbcr/fM) = 
(0.649 ± 0.004) - (2.70 ± 0.05) x 10~ 2 T/J. (b) Critical ve- 
locity normalized by the local speed of sound at the cen- 
ter of the barrier, as a function of the initial temperature. 
An exponential fit to the TWA data only gives v c /c2d = 
(1.03±0.03) exp [-(0.164 ± 0.005)T/J]. This fit is used to de- 
termine v c /c2d for the temperatures 10, 20 and 30nK, shown 
in Fig. 8. The lattice dimensions are 126 x 17 x 5 and the total 
atom number is N = 51407 atoms. The chemical potential in 
the absence of the barrier is (j,o = 2. 14 J. 



temperature dependence of the decay timescale indicates 
the importance of thermally activated processes, which 
should be measurable in experiments. 

Next, we investigate the temperature dependence of 
the critical barrier and velocity, while keeping the total 
number of atoms fixed. In Fig. 10 we plot (a) the critical 
barrier and (b) the critical velocity as a function of the 
temperature prior to phase imprinting. In addition, we 
plot the GPE prediction at T = 0. We note that the 
GPE approximation ignores quantum fluctuations, and 
is thus not the actual T = prediction. For the range 
of temperatures simulated, the critical barrier height de- 
pends approximately linearly on the temperature, as seen 
in Fig. 10(a). The gray line is a linear fit to the TWA 
data only and predicts a zero temperature critical bar- 
rier in agreement with the critical barrier height from the 
GPE simulations. 

The ratio of the critical velocity to the local sound 
speed, v c /c2D as a function of temperature is plotted on 
a log-linear scale in Fig. 10(b) and suggest an approxi- 
mately exponential relationship. Fitting only the TWA 



data to an exponential function yields a fit v c /c2D — 
(1.03 ± 0.03) exp [-(0.164 ± 0.005 )T/ J], which again is 
in agreement with the GPE result at zero temperature. 
We use this fit to estimate the critical velocity at tem- 
peratures closer to those in the experiment. We assume 
the experimental temperature to be in the range 10, 20 
or 30nK, and find that the extrapolation predicts a criti- 
cal velocity in the range u c (10nK) = (0.88 ± 0.03)c2D, 
w c (20nK) = (0.76 ± 0.03)c 2jD , w c (30nK) = (0.65 ± 
0.03)c2D- The shaded regions in Fig. 8 represent each 
of the ranges of critical velocities. This is a rough esti- 
mate for the critical velocity, based only on the data for 
= 2.14J, and does not take into account the depen- 
dence on the chemical potential which is seen in the TWA 
simulations and also weakly in the experimental data. 

We thus find that taking into account thermal fluc- 
tuations, within a TWA approach, gives a reduction of 
v c /c2D that is comparable to the one found in experi- 
ment. Further contributions to the decay could be dis- 
order in the trapping potential or other technical noise, 
that would lower v c /c2D further. 



VII. CONCLUSION 

In conclusion, we have simulated the experiment re- 
ported in [7], using a TWA simulation and, for compari- 
son, a GPE approach. We find that thermal fluctuations 
captured within TWA simulations significantly modify 
the results of GPE simulations. In particular, the crit- 
ical barrier height and the critical velocity - as defined 
in [7] - are reduced. Furthermore, by observing individ- 
ual TWA realizations, we identify the decay mechanism 
of superfluid flow in a toroidal BEC to be thermally ac- 
tivated phase slips at the barrier, which are generated 
by vortices crossing the barrier region. We also study 
the temperature dependence of the decay time scale, and 
find a strong dependence, as shown in Fig. 9. 

We compare our results with the experimental results 
reported in [7], as shown in Figs. 7 and 8. These exper- 
iments had found that v c /c e ff ~ 0.6. We find this to 
be in contradiction to GPE simulations which give ap- 
proximately v c /c e ff fa 1. However, taking into account 
thermal fluctuations within a TWA simulation offers a 
possible explanation. We find that for temperatures in 
the range 10 - 30 nK, the critical velocity can be esti- 
mated to be v c /c e ff ss 0.65-0.9. This is a reduction com- 
parable to the experimental finding. This emphasizes the 
importance of including fluctuations in the simulations of 
ultra-cold atom systems, to understand 'post-GPE' dy- 
namics. We further note that the reduction that was 
found in experiment appears to be even larger than the 
thermal reduction. This could suggest that besides the 
thermal effects that are simulated here, additional effects 
such as the visible disorder of the trap potential could re- 
duce the critical velocity to the experimentally observed 
regime. 
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Appendix A: Measuring Temperature via coupling 
and decoupling harmonic oscillators 

In order to investigate the temperature dependence of 
the superfluid decay we developed a method to measure 
the temperature of the atomic cloud in the trap. We use 
this method prior to the phase imprint, as the last step 
of the preparation stage of the numerical simulation. 

We first couple several harmonic oscillators weakly and 
adiabatically to the x- or y-component of the current, 
jx,y{?s), at different locations r s in the ring. The har- 
monic oscillator Hamiltonian takes the form: 

H ho = ^ H ° s + U hol{tn) } J Psja (r s ) 
s s 

where Hq s = \ (p 2 s +u'f lo x 2 a ) is the bare harmonic oscil- 
lator Hamiltonian and x s and p s are the position and 
momentum of oscillator s, where s — 1, . . . , A/ lo and Nh 
is the number of oscillators. 

We couple the oscillators to the current j a (r s ) = 
-iJlhr 1 (^VVs+a - V^+aVVs), where r s + a is the 
nearest-neighbor site either in x- or in y- direction. A 
schematic diagram of the harmonic oscillator thermome- 
ters is shown in Fig. 11(a). 

We turn the coupling on and off adiabati- 
cally slow, using the time dependence j{t n ) = 
{tanh [(t n - ti)/T ho \ - tanh [(*„ - t 2 )/T ho ]} /2. The 
time difference between turn on and turn off times t\ 
and ti is chosen long enough to allow the oscillators to 
equilibrate with the atomic cloud. This was checked 
by inspecting if the energy (E) = J2 S ( H 0s) /N ho had 
reached a steady state. 

The effective temperature of the cloud is calculated 
from the expectation value of the energy after decoupling, 



arctanh(?kj/j /2 (E)) 

The oscillators are initialized according to their Wigner 
distribution at finite temperature Tqj w . This corre- 
sponds to sampling from a product of Gaussian distribu- 
tions with variances a\ = 1/ [2w^ tanh(w^ o /2T ^ o )] and 
(Tp = uf lo <?x for x and p, respectively. As a further check, 
the initial temperatures of the harmonic oscillators are 
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FIG. 11. (Color online), (a) Schematic diagram of harmonic 
oscillator thermometers, (b) Effective temperature of the os- 
cillators, obtained by evaluating Eqn. (Al) at each time step. 
N = 51407 atoms. The lattice dimensions are N x = 126, 
N y — 17 and N z — 5. In this simulation 512 realizations 
are used. The turn on and off times are ti — 3150ft/ J and 
t% — 21950ft/ J; the time scale for both turning on and off the 
coupling is th = 1200ft/ J. The initialization temperatures of 
the oscillators are To,ho = 1 J, 3J and 5 J. The atomic ensem- 
ble is intialized with To/ J = 1. The harmonic trap is ramped 
on at t = 8050ft/ J, with a time constant of r = 3200ft/ J. 
After turning off the couplings to the oscillators, the temper- 
ature is estimated to be T* = (4.57 ± 0.09) J = (49.3 ± 1.0) 
nK. 



set to different values. We find that the energies of the 
oscillators converges towards the same steady state value, 
independent of the initial temperature. This is a further 
indication that the oscillators have equilibrated with the 
atomic ensembles, and that the measured temperature is 
a good estimate of the atomic ensemble temperature. 

In the numerical results reported here, six harmonic 
oscillators, equally spaced along the ring, are coupled to 
the atomic current at the trap center in the x- and y- 
directions, at the location of maximum density. Half of 
the oscillators are coupled to the x-current and the other 
half are coupled to the y-current. In Fig. 11(b) an ex- 
ample for the time evolution of the effective temperature 
of the six oscillators is plotted. For this example and 
throughout the paper, the harmonic oscillator parame- 
ters are Uho — 0.008 J, ojho = 2 J and t% = 1200ft/ J. 
These parameters were chosen in a way to minimize both 
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the computational time and the coupling strength. We 
check that the oscillators do not introduce measurable 
heating of the atoms. The green (blue) lines correspond 
to the oscillators that are coupled to the x- (y-) current. 
Each line represents an average of 512 realizations. The 
thick black line is an average over all of the oscillators 
and all of the realizations. As can be seen from the plot, 
the oscillators initialized to different temperature con- 
verge to a single temperature, the effective temperature 
stabilizes in time and the final temperature of the oscil- 
lators coupled to the x-current is indistinguishable from 
those coupled to the y-current. The oscillators initialized 
to Tq^o = 5J cool initially and then heat up again as the 
trap is ramped on. The final temperature for the data in 
Fig. 11(b) is T* = (4.57 ± 0.09) J = (49.3 ± 1.0) nK. 

We have checked that the final temperature of the atom 
cloud does not depend on the parameters chosen, within 
the error, indicating that any heating associated with 
the oscillators is negligible. We consider this method 
presented here to be generally applicable to a wide range 
of TWA simulations. 



Appendix B: Estimating the Chemical Potential 



In a homogeneous system in equilibrium and with a 
well-defined dimension, the chemical potential is a well- 
defined quantity. However, in our simulations we con- 
sider a trapped system out of equilibrium and with re- 
gions of varying dimensionality, such as the bulk and the 
barrier region. Despite this, the energy scale of a 'lo- 
cal' chemical potential is a useful quantity in discussing 
the behavior of the system, even though it only has an 
approximate meaning. 

To estimate the chemical potential for the numerical 
results, we use several approximations: (1) ^i(x), which 
estimates the chemical potential based on the density at 
location (x, 0,0), within the Thomas-Fermi approxima- 
tion; (2) the average-density- based approximation, 
which is based on the density at the trap minimum av- 
eraged around the ring; (3) /zjv, the atom- number-based 
approximation, which is calculated from the total num- 
ber of atoms, assuming a Thomas-Fermi distribution in 
both y- and z-direction; (4) /12d(x), a column-density- 
based approximation, in the limit that locally the sys- 
tem is two-dimensional (2D), calculated from the column 
density at (x, 0) and (5) ^i_d(x), based on the density in- 
tegrated over both y and z, which is applicable when 
the system is one-dimensional. The approximations ^ 
and fiN estimate the global chemical potential, while the 
others estimate the 'local' chemical potential. The ap- 
proximations fiN and fj,2D are specifically calculated to 
compare the numerical and experimental results. 

We briefly outline the method for determining each 



quantity. The energy functional of the GPE is [37]: 

h 2 



m = J 



dr 



— |VV(r)| 2 + fhMr)| 4 



(Bl) 



where V(r) is the trapping potential, V(r) — 
\m{yj^y 2 + u> 2 z 2 ). To determine the equilibrium state 
the GPE energy functional is minimized, which gives 

[-h 2 /{2m)\V^{v)\ 2 + ,g|^(r)| 2 + V(r)] ^(r) = ^(r), 

where fi is introduced as a Lagrange multiplier. Within 
the Thomas-Fermi approximation, the kinetic energy 
term is neglected and the density is 

n(r) = |^(r)| 2 = 5 - 1 ( M -y(r)). (B2) 

At the trap minimum, the local chemical potential is 

MO) = gn(x, 0,0). 

This approximation applies when the chemical potential 
is greater than the trapping energies, jj, > ftio yi hcu z . 

To improve this estimate in the numerical evaluation, 
we calculate the average density at the trap minimum, 
where V(r) = 0, 

X 

We calculate this quantity in the absence of the barrier, 
to estimate the bulk chemical potential. 

The chemical potential can also be determined by sum- 
ming the density over all space in Eqn. (B2), using the 
condition J drn(r) — N, and solving for the chemical po- 
tential as a function of the total number of atoms. The 
total-number-based approximation for the chemical po- 
tential is 



{ gNmu y uj 2 
Mat = I — 



TtL„ 



1/2 



^ UlNh 2 uj y uj z 



2irL T J 



1/2 



This expression assumes a Thomas-Fermi profile in y- 
and z-direction. This quantity was used in the experi- 
ment in Ref. [7], and is used here to compare the simu- 
lations and the experimental results in Fig. 7. 

Furthermore, the chemical potential at the barrier was 
calculated from the two-dimensional column density at 
the barrier in Ref. [7]. At the barrier, the density is suffi- 
ciently suppressed that the local chemical potential is less 
than the harmonic confinement in z-direction, so that it 
is effectively 2D. The wavefunction is separable, ip(r) = 
tpTF(x,y)ipho(z), into a product of the harmonic oscil- 
lator ground state tpho(z) = (l/irl 2 ) 1 ^ exp (— z 2 /2l^), 
where l z = yfh/mu) z , and the Thomas-Fermi distri- 
bution in x- and y-direction. The harmonic oscillator 
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ground state is normalized as J dz\tpho(z)\ 2 = 1 so that 
\iPtf(x, y)\ 2 corresponds to the column density measured 
in the experiment. We substitute ipTF(x,y),ipho{z) into 
Eqn. (Bl) and integrate over z to get 



in which several limits for the chemical potential were 
discussed. 

For a Bose condensate in a cylindrical geometry, the 
phonon velocity is approximately given by 



m = J 



dxdy 



+ 



2m 
9 1 



2 y/2nl 



-\lpTF\ 



Again, we ignore the kinetic energy term, and subtract 
the constant offset to the chemical potential due to the 
harmonic oscillator energy, fuv z /2. We minimize the total 
energy, while fi controls the total number of particles. 
The resulting 2D column density is 



n2D(x,y) = \^tf\ 2 = — 

92D 



where g 2 u = g/V2-rrl z . Solving this expression for the 
chemical potential in terms of the peak column density, 
U2D, which occurs at y = 0, yields 

H2d(x) = g2Dn 2 D(x,0). 

The chemical potential for the quasi- ID case can be de- 
termined in a similar manner, by replacing the full wave- 
function with V( r ) — *P{ x )ipho,y(y)ipho,z(z)- The result- 
ing chemical potential is 

Hid(x) = giDn 1D (x) 

where gm — g/{2irl z l y ), and nm is the density inte- 
grated along y- and z-direction, ni£>(x) = J2 y z n i x , V> z )- 



Appendix C: Speed of Sound 

An important dynamic quantity of a condensed Bose 
gas is the phonon velocity. For a homogeneous, weakly 
interacting Bose gas in 3D it is given by c s = \J gnjm. 
However, the system that we consider here has a spatially 
varying density, even to the degree that the dimension of 
the system varies, for example in the vicinity of the bar- 
rier. We therefore introduce several limiting expressions 
for the phonon velocity similar to the previous section, 



C3d(x) 



gn(x, 0, 0) Jl 



2m 



j h 



The density n has been replaced by the average density 
over the cylinder, n = n(x, 0,0)/2, where we assume a 
Thomas-Fermi profile in y- and ^-direction, see e.g. [30, 
38]. This is valid when the chemical potential is greater 
than the harmonic confinement energy in the transverse 
directions. When hu) y < ^i(x) < ftw Zl the system is quasi- 
2D and the local speed of sound is 



c 2 d(x) = 



2g2Dn 2 D(x,0) Ah2d{x) Jl 



3m 



3J 



This can be derived from the low-energy excitation spec- 
trum within the Bogoliubov de Genes approximation, 
as in Refs. [38, 39], starting from the 2D Hamiltonian. 
It was obtained in Ref. [31] within a hydrodynamic ap- 
proach. The key step is to take the average over the 
column density, n 2D (x,y) = (^ 2 d - \muly 2 ) / 'g 2D over 
y- n 2D = J dyn 2D (x,y) = 2n 2D (x,0)/3. 

At the critical barrier, the condition ftw y < fj,(x) < hu z 
is fulfilled for all of the simulations using the truncated 
Wigner approximation reported in this paper. Thus we 
use c 2 d as the best approximation at the barrier. Addi- 
tionally, the experimental data is rescaled and the critical 
velocity is compared with c 2 ]j instead of c 3 £> , as was orig- 
inally done in [7]. 

When fj,(x) < huj y , haj z , the system is quasi-lD and the 
local speed of sound is the same as in the homogenous 
case, except with the ID interaction parameter and ID 
density, 



cid(x) = 



g\Dn 1D {x) _ 2m D {x) Jl 



m 



J 



The system is quasi- ID at the critical barrier for some 
of the GPE runs, in particular for the case of = 2.2 J. 
We then compare the local velocity to c\d for the GPE 
data Fig. 2. 

We note that the cross-over from the quasi-lD to the 
quasi-2D regime is not a sharp transition. We find that 
in the intermediate regime of fj,(x) w Huj y , cm is typically 
<~ 10-15% lower than c 2 _d. 
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